% Capillary Blood Flow Analysis 
%  based on TK's adoption of BC/VN's code. This is written to work with
%  drawRBCKymline_batch and batchKymograph to do all this processing in
%  matlab along with drift correction

% Set Parameters Based on Image Acquisition Settings
%  Inputs:
%  inputimg = linescan/space-time image with space in X-axis and time in Y-axis
%  showimg = display image comparing the effect of elementwise demeaning, vertical demeaning and 3x3 vertical Sobel filtering
%  saveimg = save image in various formats (.fig,.jpg,.eps,.ai)
%  delx = DeltaX, microns/pixel
%  delt = DeltaT, ms/line
%  hi = height of image segment to be processed in pixels
%  lineskip = number of lines before next image segment starts
%  xrange = 2-element vector specifying range of pixels in space-dimension to use in the image segment
%
% updates from v2.0:
% set up for the visual stim paradigm - a bit hacky right now but can
% automate more in the future
function results = capillaryBloodFlow_analysis_v2_1(filenames,fighandle,objective)

if nargin > 1 && ~isempty(fighandle)
    doplot = true;
else
    doplot = false;
end
if nargin < 3 || isempty(objective)
    objective = '25x';
end

metaname = strrep(filenames(1).name,'_kymograph.mat','.tif');       %assumes the same imaging parameters for all the kymographs; can update this in future iterations
meta = getSImetadata(metaname);

delx = pixel2micron_2p(meta.zoom,objective);
delt = 1000/meta.framerate; %ms/linescan - note this is NOT how long it takes to aquite one line of the x-line tall image;

%parameters to set the windowing of the kymograph (see hybridvel_sobelonly)
hi = 100; 
lineskip = 25;  

%% load in kymographs generated from batchKymograph.m
numkymos = length(filenames);
kymographs = cell(size(filenames));
for n = 1:numkymos
    load(filenames(n).name,'kymo');
    kymographs{n} = kymo;
    clear kymo
end

%%  Sobel Filtering & Iterative Radon Transformation - Find Velocity Over Time
%  Outputs:
%  angle: [angle(deg), minstep(deg), location(pixel), dels1(deg), deln(deg), %dv/v, iter, irl, speed (mm/s)]
%  utheta: angles used to process the image segments
%  uvar: variance on Radon transform at each angle at each image segment

time = [];
velocity = [];
baseline_velocity = [];
percent_velocity_change = [];
velocity_sem = [];

% the results time conversion is set by the framerate of the system + the
% lineskip parameter. 
baselineStimTime = [3,5];   %in seconds, the range of time used for calculating the pre-stimulation baseline
baselineframes = unique(round((round(baselineStimTime(1)*meta.framerate):round(baselineStimTime(2)*meta.framerate))/lineskip));
for i=1:numkymos
    
    currKymo = kymographs{i}*100;
    imsize = size(currKymo);
    xrange = [1 imsize(2)];
    
    % Hybrid Velocity Function Call
    warning off; % to suppress... "Warning: The EraseMode property is no longer supported and will error in a future release."
    [angle,utheta,uvar] = hybridvel(currKymo,[],[],delx,delt,hi,lineskip,xrange);
    warning on;
    
    % Smoothing, Calculate Percent Change In Velocity Relative To Baseline
    time(:,i) = angle(:,3,3)*delt/1000;
    velocity(:,i) = smooth(angle(:,9,3), 15);
    baseline_velocity(:,i) = velocity(baselineframes,i);%these indices correspond to 2 seconds before stim, which is used in Chow & Nunez et al 2020 for normalization
    avg_baseline = mean(baseline_velocity(:,i));
    percent_velocity_change(:,i) = ((velocity(:,i)-avg_baseline)/avg_baseline)*100;
     
end

results.time = time(:,1);
results.velocity = velocity;
results.baseline_velocity = baseline_velocity;
results.percent_velocity_change = percent_velocity_change;

if doplot
    figure(fighandle);
    cla;
    plot(time(:,1), mean(percent_velocity_change,2), ['-o' 'b' ]);%plot average of kymographs
    xlabel('time (s)'); ylabel('% delta V / V baseline'); 
    yl = get(gca,'ylim');
    hold on
    h = fill([5,15,15,5],[yl(1),yl(1),yl(2),yl(2)],[.5,.5,.5]);
    h.FaceAlpha = .3;
    h.LineStyle = 'none';
    hold off
    pause(.1);
end



